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Hydrodynamics merged with single-freeze-out statistical hadronization is used to describe the 
midrapidity hadron production in relativistic heavy-ion collisions at the highest RHIC energies 
{^/snn = 200 GeV) and to make predictions for the LHC (y^sjvjv = 5.5 TeV). Thermodynamic 
properties of the high-temperature strongly-interacting quark-gluon plasma are taken from lattice 
simulations, at low temperatures the hadron-gas equation of state is used, while in the cross-over 
region an interpolation between the two equations of state is constructed. Boost invariance is 
assumed for the midrapidity calculations. The initial condition for hydrodynamics is obtained from 
a Glauber profile for the entropy, with the initial central temperature Ti. The conditions obtained 
from the hydrodynamic expansion at the freeze-out temperature T/ are used as input for the thermal 
event generator THERMINATOR, which accounts for a complete treatment of hadronic resonances. Basic 
physical observables are obtained: the transverse-momentum spectra, the elliptic flow coefficient V2, 
and the HBT radii. The femtoscopic observables are evaluated with the help of the two-particle 
method which accounts for the resonance decays and Coulomb final-state interactions. The problem 
of a simultaneous description of all discussed observables is addressed, with the conclusion that 
at the highest RHIC energies our approach gives a quite satisfactory global description of soft 
hadronic observables, which agree with the data at the level of 10-15%. Some discrepancies may be 
attributed to the absence of the final-state elastic interactions among hadrons. Extrapolating Ti to 
higher values allows for global predictions for soft hadronic physics at the LHC. We test Ti = 400, 
450, and 500 MeV, and observe the expected growth of particle multiplicities and the increase of 
the flow, resulting in smaller slopes of the pr-spectra. The elliptic flow of pions exhibits saturation, 
with V2 remaining practically constant, while the HBT radii increase moderately with Ti. 
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I. INTRODUCTION 

Hydrodynamics has become the standard framework 
for the description of the intermediate stages of relativis- 
tic heavy-ion collisions [11,0,0,0,0]. In this paper, en- 
couraged by the success of this approach in describing 
some of the basic RHIC data @, 0, & 0, E, EI El , 
we use the recently developed hydrodynamic approach 
01,0 merged with the statistical-hadronization 
Monte-Carlo model to globally describe the midra- 
pidity hadron production at the highest RHIC energies 
and to make predictions for the future experiments at 
the Large Hadron Collider (LHC). Our approach uses 
standard methods of the field, however we take an effort 
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to combine the best features for all entering elements of 
the existing analyses, except for the absence of elastic 
rescattering in the final state. We attempt a global fit 
to all main midrapidity observables, i.e., the transverse- 
momentum spectra, the elliptic-flow coefficient, as well 
as the pionic HBT radii. The main outcome is that for 
pions and kaons a global fit works at the level of 10-15%, 
which in our view is satisfactory, baring in mind system- 
atic uncertainties in various elements of the approach, 
such as the determination of the initial condition, the 
lack of detailed knowledge of the equation of state, the 
assumed boost invariance, etc. The quality of the global 
fit is calling for predictions for the LHC energies, which is 
the main topic of this work. The predictions of other hy- 
drodynamic models for phenomena expected at the LHC 
have been recently summarized in Ref. [Isj ]. 

Our implementation of inviscid hydrodynamics 
(Sect. Ull) incorporates the state-of-the-art knowledge of 
the hadronic equation of state. The relevant quantity 
is the sound velocity, c s , considered as a function of 
temperature. The value of c s for the high-temperature 
strongly-interacting quark-gluon plasma is taken from 
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the recent lattice simulations [ijj, which exhibit a sub- 
stantial departure from the ideal-gas Stefan-Boltzmann 
limit even at temperatures significantly above the 
cross-over temperature T c ~ 170 MeV. In particular, 
c 2 s is significantly below the ideal-gas value of 1/3. At 
low temperatures the hadron-gas equation of state is 
used, with a complete treatment of resonances [13] ■ In 
the limit T — > the function c 2 s {T) approaches zero 
as (r/mTr) 1 / 2 , characteristic of the massive pion gas 
at low temperatures. In the cross-over region near 
T c a simple-minded interpolation between the high- 
and low-temperature equations of state is constructed. 
In accordance to the present knowledge, no phase 
transition but a smooth cross-over is built in, thus c s 
does not drop to zero at T c . With the assumed features 
our hydrodynamics leads to smooth, laminar flow; no 
shock-waves are formed. 

Boost invariance is assumed for the considered midra- 
pidity analysis, which is a good approximation for the 
highest RHIC and LHC energies at the rapidity window 
\y\ < 1 [21L l22l|. This assumption reduces the number of 
independent space dimensions, which greatly simplifies 
the numerical analysis. The hydrodynamic equations are 
solved with the technique described in Refs. [3, EH, EH] 
which is a generalization of the method introduced by 
Baymetal. in Ref. [23]. The initial condition for the evo- 
lution is obtained from a Glauber profile (Sect. [nil, with 
the initial central temperature Tj. The evolution pro- 
ceeds until the temperature drops down to the freeze-out 
temperature Tf. These two temperatures are the free pa- 
rameters of the approach. Entropy conservation is used 
as a numerical test, satisfied at the relative level of 10~ 4 . 
The freeze-out hypersurface and the flow profile obtained 
from the hydrodynamic expansion are then used as input 
for the thermal event generator THERMINATOR [t?| . This 
code implements the statistical hadronization, account- 
ing for a complete treatment of hadronic resonances (the 
included resonances and their branching ratios are the 
same as in the SHARE package [2J]). Rescattering after 
the chemical freeze-out is not incorporated, which is an 
approximation working fairly well for the pions, as may 
be inferred from the results presented in Ref. [ijj]. We 
note that a Monte Carlo generator of functionality simi- 
lar to THERMINATOR has been released recently (251 ]. 

Basic physical observables are calculated from the sam- 
ple of events generated by THERMINATOR (Sects. El EI 
and IVIip . For the highest RHIC energies the pion and 
kaon transverse-momentum spectra agree well with the 
data, with a similar quality matching of the elliptic flow 
coefficient Vsipr)- For the protons the model spectra are 
somewhat too steep and V2 too large, which is proba- 
bly due to the absence of the final-state elastic rescat- 
tering. The femtoscopic observables are evaluated with 
the help of the two-particle method accounting for the 
effects of resonance decays and the Coulomb final-state 
interactions (Sects. [Vl IVII and lVIIj) . The use of the two- 
particle method imitating the experimental analysis is a 
clear advantage over the use of simple parameterizations 



of the emitting source (26J. The RHIC pionic HBT radii 
are reasonably reproduced in our approach, with i? s id e 
10-15% below the experimental data, i? ou t within the ex- 
perimental errors, and i?i on g about 15% too large. The 
ratio i?out/-Rsidc is about 1.2-1.25 and almost constant as 
a function of the the transverse momentum of the pair, 
which is still away from the experimental ratio, but con- 
siderably better than in many other hydrodynamic ap- 
proaches. All in all, at the highest RHIC energies the 
approach gives in our view a quite satisfactory global de- 
scription of the soft pionic observables, which agree with 
the data at the level of 10-15%. Some discrepancies may 
be attributed to the final-state elastic interactions among 
hadrons, not included in our description. 

Extrapolating Tj to higher energies allows us to make 
global prediction for the soft hadronic physics at the 
LHC, which are presented in Sect. IVIll We start the hy- 
drodynamic evolution from a higher initial temperatures, 
using a few values in the expected range: T% = 400, 450, 
and 500 MeV, and keep all other parameters fixed, in par- 
ticular Tf — 150 MeV. We observe the expected growth 
of particle multiplicities and the increase of the flow, re- 
sulting in smaller slope of the pr-spectra. The elliptic 
flow of pions does not increase any more, showing the ex- 
pected saturation of v%, and even decreases slightly at the 
highest tested temperature of T, = 500 MeV. The HBT 
radii increase moderately with the temperature: i? s ide 
and -Riong roughly 1 fm for each 100 MeV increase of T,-, 
and i? ou t even slower. We provide some simple formulas 
describing the change of the considered quantities with 
Ti. 

Finally, we remark that THERMINATOR may be used 
straightforwardly to test the soft physics in detector sim- 
ulations at LHC. It produces particles with full decay 
history which can be fed directly to transport codes. It 
includes non-trivial and physically well-motivated pre- 
dictions for flow phenomena (both the radial and ellip- 
tic flow) with particle type dependence naturally built 
in. It also provides information on emission points for 
femtoscopic simulation. The non-trivial, dynamic (con- 
taining space-momentum correlations) emission function 
allows for the testing more advances femtoscopic analy- 
sis techniques, e.g. spherical harmonics decompositions 
or imaging [23, HE HE H|- 



II. EQUATION OF STATE AND 
PERFECT-FLUID HYDRODYNAMICS 

In our studies we focus on the description of the midra- 
pidity region. Statistical models applied to the highest- 
energy RHIC data yield the baryon chemical potential 
HB ~ 20 - 30 MeV [Hl|, \M, M, 0, El, while the predic- 
tions for LHC energies give [Ib ~ 0.8 MeV [HH]. Hence, 
for the hydrodynamic equations we can approximately 
assume that the baryon chemical potential vanishes. In 
this situation the whole information on the equation of 
state is encoded in the temperature-dependent sound ve- 
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FIG. 1: (Color online) The square of the sound velocity as a 
function of temperature. The high-temperature part comes 
from the lattice simulations of Ref. (jg| [. the low-temperature 
part from the hadronic gas with resonances, while the tran- 
sition region near the cross-over is obtained with simple in- 
terpolation. Here we take T c = 170 MeV as the scale for 
T. 




T [MeV] £ [GeV / fm 3 ] 



locity c s (T) In the hydrodynamic calculations we 

use the function c s (T) introduced in Ref. and la- 
beled there as "case I". At low temperatures it is given by 
the sound velocity of a hadron gas (with a complete set 
of hadronic resonances) . The function c 2 (T) approaches 
zero as (T/m^) 1 ^ 2 , characteristic of the pion gas. At 
high temperatures the equation of state coincides with 
the recent lattice simulations [lj|. In the intermediate- 
temperature region we take the simple smooth interpola- 
tion between the hadron gas and lattice results, see Fig.[TJ 
As mentioned above, the knowledge of the function c s (T) 
allows us to determine all other thermodynamic proper- 
ties of our system. In Fig. [2] we display the corresponding 
entropy and energy densities as functions of T, as well as 
the pressure and the sound velocity as functions of the 
energy density. The region in the energy density where c 2 
drops to zero is exceedingly small and invisible in Fig. [2j 
while the full dependence can be seen in Fig. [TJ Other 
equations of state and their impact on the physical ob- 
servables were recently studied in Ref. (37l |. 

We stress that by using the lattice results we take into 
account the non-perturbative aspects of the plasma be- 
havior, which may be regarded as the effective inclusion 
of the strongly-interacting quark-gluon plasma; large de- 
viations from the ideal-gas behavior directly indicate the 
non-negligible interactions present in the plasma. In par- 
ticular, c 2 is significantly below the ideal-gas value of 1/3 
also at temperatures way above T c . Note that in agree- 
ment with the present knowledge, no phase transition is 
present in the system, but a smooth cross-over, there- 
fore c s does not drop to zero at T c but remains a smooth 
function. Also, the criterion for absence of shock-waves 
is satisfied [H, [38|, which makes the evolution simple to 
solve numerically - we use the adaptive method of lines in 
the way implemented in the MATHEMATICA package. 
The space directions are discretized, and the integration 
in time is treated as solving of the system of ordinary 



FIG. 2: (Color online) The two left panels: the volume densi- 
ties of entropy and energy, scaled by T 3 and T 4 , respectively, 
shown as functions of the temperature. The two right pan- 
els: the pressure and sound velocity shown as functions of the 
energy density. The presented thermodynamic functions fol- 
low directly from the temperature-dependent sound velocity 
shown in Fig. [T] 



differential equations with the adaptive time step. En- 
tropy conservation is used as a numerical test, which is 
satisfied at the relative level of 1CP 4 . 

With the assumed features our hydrodynamics leads 
to smooth, laminar flow. In Ref. [16j it was argued that 
the use of the sound- velocity function depicted in Fig. Q] 
is physically attractive, since it leads to relatively short 
time scales of the hydrodynamic evolution. This, in turn, 
helps to obtain the satisfactory description of the HBT 
radii at the highest RHIC energies. 

The relativistic hydrodynamic equations of the perfect 
fluid follow from the energy-momentum conservation and 
the assumption of local equilibrium. For the baryon-free 
matter they read 

u"0 M (2V) = d"T, (1) 
a^au") = 0, (2) 

where T is the temperature, s the entropy density, and 
u^ 1 = 7(1, v) denotes the four-velocity of the fluid. The 
normalization condition u fJ u /1 = 1 implies that only three 
out of four equations appearing in |(T]) are independent. 
As mentioned above, in the analysis of the evolution in 
the midrapidity region we additionally assume that the 
system is boost-invariant. This assumption introduces 
another constraint, hence, in this case Eqs. {J) and ^ 
reduce to three independent equations, which may be 
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FIG. 3: (Color online) The notation for the position and ve- 
locity of the fluid element in the transverse plane. 



written in the following form [3£ 
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dt 
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dr 
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(rTjv) 



(rts'yv cos a) + 
d 

r cos a— (Tj) 
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sma 



dT 
dr 



— (tsjvsma) 
9 

sinaA (T7) 
cos a dT 



r 



(3) 



Here t, r=\Jx 2 



y 2 , and <f) = tan _1 (y/x) are the time 
and space coordinates which parameterize the transverse 
plane z = (for boost-invariant systems, the values of 
all physical quantities at z ^ may be obtained by the 
Lorentz transformation) . The quantity v is the transverse 
flow, 7 = l/y/l — v 2 is the Lorentz factor, and a is the 
dynamically determined angle between the direction of 
the transverse flow and the radial direction, see Fig. 03 
The differential operator d/dt represents the complete 
time derivative and is defined by the formula 



d d d v sin a d 

dt dt dr r dS 



(4) 



Equations (E3) are three equations for four unknown func- 
tions: T, s, v, and a. To obtain a closed system of equa- 
tions one needs to supplement them with an equation of 
state, i.e., with a relation connecting T and s. Alter- 
natively, one may fix the temperature-dependent sound 
velocity, 



c 2 s (T) 



dP 

~d7 



s dT 
T~ds' 



(5) 



to a known function. This method is advantageous, as 
it directly uses the available information on the sound 
velocity. 



III. INITIAL CONDITIONS 

Similarly to other hydrodynamic calculations, we as- 
sume that the initial entropy density of the particles pro- 
duced at the transverse position point xt is proportional 



to a profile obtained from the Glauber approach. Specifi- 
cally, we use a mixed model [13, El| , with a linear combi- 
nation of the wounded-nucleon density pw ( x t) and the 
density of binary collisions pbin (^r), namely 



s (xt) oc p (xt) 



I - 



■ pw (Xt) + Kp hin (xt) ■ (6) 



The case k = corresponds to the standard wounded- 
nucleon model [Hj, while k — 1 would include the binary 
collisions only. The PHOBOS analysis [4l| of the particle 
multiplicities yields k = 0.12 at ^/snn = 17 GeV and 
k = 0.14 at ytsNN — 17 GeV. In this paper we assume 



0.14 



(7) 



for both the highest RHIC and the LHC energies. Since 
the density profile from the binary collisions is steeper 
than from the wounded nucleons, increased values of 
k yield steeper density profiles, which in turn result in 
steeper temperature profiles. 

In the hydrodynamic code, the initial conditions are 
specified for the temperature profile which, according to 
Eq. jBJ, takes the form 



T(n n it,x T ) = T s 



p(xt) 
P(0) 



(8) 



where T$(s) is the inverse function to the function s(T), 
and Si is the initial entropy at the center of the sys- 
tem. The initial central temperature Ti equals Ts(si). 
Throughout this paper we take the initial time for the 
hydrodynamic evolution to be 



Tinit = 1 fm. 



(9) 



The wounded-nucleon and the binary-collisions densi- 
ties in Eq. ([8]) are obtained from the optical limit of the 
Glauber model, which is a very good approximation for 
not too peripheral collisions [43] . The standard formulas 
are El 



pw (xt 




and 



Pbin (Xt) = oT A 



x T T A 



b 

2 



X T 



(11) 



In Eqs. (fTfjj) and (fir) b is the impact parameter, a 
is the nucleon-nucleon total inelastic cross section, and 
Ta (x, y) denotes the nucleus thickness function 



T A (x,y) = J dz 



p(x,y,z) . 



(12) 
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FIG. 4: (Color online) The parameterization of the freeze-out 
hypersurface. 

For RHIC energies we use the value a — 42 mb, while for 
LHC we take a — 66 mb. The function p(r) in Eq. lfl2|) 
is the nuclear density profile given by the Woods-Saxon 
function with the conventional choice of parameters: 

po = 0.17 far 3 , 
r Q = (1.12A 1/3 -0.86A- 1/3 )fm, 
a = 0.54 fm. (13) 

The atomic number A is 197 for RHIC (gold nuclei) and 
208 for LHC (lead nuclei) . The value of the impact pa- 
rameter in Eqs. (fit))) and (flTTl depends on the considered 
centrality class. 

We stress that the shape of the initial condition JH]) 
is important, as it determines the development of the 
radial and elliptic flow, thus affecting such observables 
as the £>T-spectra, v<i, and the femtoscopic features. On 
qualitative grounds, sharper profiles lead to more rapid 
expansion. Several effects should be considered here. 
Firstly, as discussed in Ref. jiij], hydrodynamics starts 
a bit later, when the profile is less eccentric than origi- 
nally due to early evolution of partons in the pre-hydro 
phase. On the other hand, statistical fluctuations of 
the axes of the second harmonic in the distribution of 
Glauber sources (wounded nucleons, binary collisions) 
[S EE EE ES EE S3, ED, E3 leads to a significant en- 
hancement of the eccentricity, especially at low values 
of the impact parameter. Thus the initial eccentricity 
may in fact be smaller or larger than what follows from 
the application of the Glauber model. This contributes 
to the systematic model uncertainty at the level of, say, 
10-20%. This uncertainty could only be reduced by the 
employment of a realistic model of the pre-hydrodynamic 
evolution and is outside of the scope of this work. With 
this uncertainty in place, one should not expect or de- 
mand a better agreement with the physical observables 
than at the corresponding level of 10-20%. 



IV. STATISTICAL HADRONIZATION WITH 
THERMINATOR 

In order to calculate the physical observables we first 
use the hydrodynamic code with the described initial con- 
ditions to determine the freeze-out hypersurface S. In 
our approach it is defined by the condition of a constant 
freeze-out temperature Tf, and parameterized in the gen- 
eral way by the following equations 

t = d (<fi, C, 77) sin£ cosh 77, z = d (4>, C, if) sin £ sinh 77, 
x — d (</), (, 77) cos £ cos (p, y = d (0, £, 77) cos £ sin <j>. 

(14) 

The variable 77 is the space-time rapidity defined as usual 
by the formula 

77=-ln^— =tanh-(-). (15) 

The parameterization (fT4| yields the compact expres- 
sions for the proper time r and the transverse distance 
r, 



t = \/t 2 — z 2 = d ((f), C, 77) sin C, 

r = \/ x 2 + y 2 = d (4>, C, 77) cos£- (16) 

At any given value of the space-time rapidity 77 the po- 
sition of the point on the hypersurface S is defined by 
the two angles, <j) and C, and the distance from the origin 
of the coordinate system, d((j>,C,r)). The angle <fi is the 
standard azimuthal angle in the y — x plane as used in 
the hydrodynamic equations ([3]), while the angle £ is the 
azimuthal angle in the t — r plane. We have introduced 
the angle ( because in most cases the freeze-out curves in 
the r — r plane may be treated as functions of this param- 
eter. The use of the transverse distance r is inconvenient, 
since very often two freeze-out points correspond to one 
value of r, see Fig. |H 

With the standard definition of the four-momentum in 
terms of the rapidity and transverse momentum, 

j/ 1 = (rriT cosh y,pr cos 4> p -,pt sin <p p , mj- sinh y ) . (17) 



where mx = y/m 2 + p 2 -, is the transverse mass, and with 
the standard definition of the element of the hypersur- 
face dS M , we find the explicit form of the Cooper- Frye 
integration measure [53] 
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dE^p^ = d 2 sm( 



dcos( (rriT sin ( cosh (77 — y) + p T cos £ cos (</> — </> p )) 



+ ^ cosC (-wit cos C cosh (77 - y) + pTsin£cos ((f) - 4> p )) + ^prsin ((f) - P ) +^ cot (mTsinh (77 - ?/) 



(18) 
drjdtfidC,. 



This equation, when used in the Cooper- Frye formula [531], leads to the six-dimensional particle density at freeze-out 



dN 



dyd(j) p pT dpT di]d<j)d( 



d 2 sin ( 
' (2tt) 3 



dcos £ (tut sin £ cosh (i] — y) + pr cos £ cos ((f) — <fi p )) 



dd „ , t , \ . > / , , \s dd 

+ — cos C (-771T cos C cosh (7]- y) + pr sm 4 cos (0 - P J) + — pr sin ( 



x < exp 



: cosh(y — 7/) ji = cos(0 + a — P ) — p/i 



± 1 



(19) 



r 



In the transition from (fT8| to (fT9| we made use of the 
boost-invariance and removed the term dd/dr/, since the 
function d depends only on (f> and C- The last line in lfl9|) 
contains the equilibrium statistical distribution function, 
with —1 for bosons and +1 for fermions. The argument 
is the Lorentz-invariant product p^u p . The parameter 
/3 = 1/T is the inverse temperature and g denotes the 
spin degeneracy factor. The four-velocity field has been 
expressed in terms of the transverse flow v and the dy- 
namical angle a which depend on the space-time posi- 
tions on the freeze-out hypersurface. The statistical dis- 
tribution functions include the chemical potential, 



/j = iibB + ll s S + /J/3/3 



(20) 



with -B, S 7 and ^3 denoting the baryon number, 
strangeness, and isospin of the particle, respectively, 
while the baryonic, strange, and isospin chemical poten- 
tials assume the values: /is = 28.5 MeV, fis = 9 MeV, 
//j 3 = -0.9 MeV at the highest RHIC energies 0, and 
Ms = 0.8 MeV and fj, s = Hi 3 = MeV at the LH C 0. 

In order to describe the statistical hadronization and 
production of particles we first run our hydrodynamic 
code and extract the functions describing the freeze-out 
hypersurface and flow: d — d(<f>,Q, v = v((f), £), and 
a = a(cj), Q. In the next step these functions are used as 
input for the THERMINATORcode, which generates, accord- 
ing to the formula l|19p , the distributions of the primor- 
dial particles. The primordial particles include the stable 
hadrons as well as all hadronic resonances. The reso- 
nances decay through strong, electromagnetic, or weak 
interactions at the random proper time controlled by the 
particle's life-time, and at the location following from the 
kinematics. As the final result we obtain the stable parti- 
cle distributions including the feeding from the resonance 
decays. Since the memory on all decays is kept, one may 
apply the experimental cuts or the weak-decay feeding 
policy, which facilitates the more accurate comparison of 
the model to the data. For example, one may extract 
easily the model proton spectrum which does not include 



the feeding from the A decays. 

The momentum distribution functions are obtained by 
integrating over </>, £, and 77 



dN 
dyd 2 pq 



7T 27T 



= / d V / dC 



dN 



dyd(j) p pT dpx drjdtfid^ 



(21) 



For cylindrically asymmetric collisions and midrapidity, 
y = 0, the transverse-momentum spectrum has the fol- 
lowing expansion in the azimuthal angle of the emitted 
particles 



dN 



dN 



dyd 2 pT dy 2irpT dpx 



(1 + 2v 2 {pt) cos(2<P P ) + ...). 



(22) 

Eq. (|22| defines the elliptic flow coefficient v 2 . 

The well-known problem of the freeze-out condition 
in the Cooper-Frye formulation is that the hypersur- 
face from hydrodynamics typically contains a non-causal 
piece, where particles are emitted back to the hydrody- 
namic region In our approach we follow the usual 
strategy of including only that part of the hypersurface 
where dE^p^ > 0. To estimate the effect from the non- 
causal part we compute the ratio of particles flowing 
backwards (i.e. where < 0) to all particles. For 
the cases studied in the following sections we find this 
ratio to be a fraction of a percent. Such very small val- 
ues show that the known conceptual problem is not of 
practical importance for our study. 



V. PION CORRELATION FUNCTION 

The correlation function for identical pions is obtained 
with the two-particle Monte-Carlo method discussed in 
detail in Ref. [H, [HfiJ . In this approach the evaluation of 
the correlation function is reduced to the calculation of 
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the following expression 
C(q,k) = 

E E $A(q - Pi + Pj )5 A (k - ±(p t+Pl ))\y(k*,r*)\ 2 



E E 5A(q-Pi+Pj)5A(k - ?(Pi+Pj)) 



where 6a denotes the box function 
Sa(p) 



i if W<#,kl<#,b*|<# 







2 > \Fy\ — 2 ' — 2 

otherwise. 



(23) 



(24) 



In the numerator of Eq. f23j) we include the sum of the 
squares of modules of the wave function calculated for all 
pion pairs with the relative momentum q (we use the bin 
resolution A = 5 MeV) and the pair average momentum 
k. For non-central collisions we only provide azimuthally 
integrated HBT radii. For each pair the wave function 
V?(fe*,r*), including the Coulomb interaction, is calcu- 
lated in the rest frame of the pair; k* and r* denote the 
relative momentum and the relative distance in the pair 
rest frame, respectively. In the denominator of Eq. lj23j) 
we put the number of pairs with the relative momentum q 
and the average momentum k. The correlation function 
(|23| is then expressed with the help of the Bertsch-Pratt 
coordinates hp, q ut, 9side, (/long and approximated by the 
Bowler-Sinyukov formula [57l.l58| 



C(q,k) = (l 

_ r>2 2 
-'''side "side 



A) + \K cou i(q inv ) [l + exp (-Rl 



l/io 



)] 



(25) 



where K cou i(q mv ) with q mv = 2k* is the squared Coulomb 
wave function integrated over a static gaussian source. 
We use, following the STAR procedure [59], the static 
gaussian source characterized by the widths of 5 fm in 
all three directions. Four kx bins, (0.15 — 0.25), (0.25 — 
0.35), (0.35-0.45), and (0.45-0.60) GeV, are considered. 
The 3-dimensional correlation function with the exact 
treatment of the Coulomb interaction is then fitted with 
this approximate formula and the HBT radii i? ut, -Rsidc, 
and i?iong are obtained. They can be compared directly 
to the experimental radii. 
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FIG. 5: (Color online) In-plane and out-of-plane freeze-out 
curves, i.e., the intersections of the freeze-out hypersurface 
with the planes y = and x = 0, obtained for central RHIC 
collisions; 6 = 1 fm, T t = 320 MeV, and T f = 150 MeV. The 
two curves overlap, indicating that the system at freeze-out 
is almost azimuthally symmetric in the transverse plane. 




FIG. 6: (Color online) Transverse-momentum spectra of 7r + , 
K + , and protons. The PHENIX experimental results [60] for 
Au+Au collisions at ^/snn = 200 GeV and the centrality 
class 0-5% (points) are compared to the model calculations 
(solid lines) with the same parameters as in Fig. 



we take (for simplicity) the same value of the impact pa- 
rameter, b = 7 fm 1 . 



VI. DESCRIPTION OF THE 
HIGHEST-ENERGY RHIC DATA 

In this Section we apply our model to describe the soft 
hadron production in the relativistic heavy-ion collisions 
studied at the highest RHIC energy, ^/snn = 200 GeV, 
by the PHENIX [6(1 [HI and STAR Collaborations [H§|. 
We consider the most central events, given by the cen- 
trality class c = - 5%, where we use the impact param- 
eters 6=1 fm. We also use the data from the centrality 
classes c = 20 - 30% (transverse-momentum spectra and 
the HBT radii) and c = 20 - 40% (elliptic flow), for which 



A. Central collisions 

First, we consider the centrality class - 5% with the 
corresponding value of the impact parameter 6=1 fm. 



1 In a detailed analysis a more accurate study in each centrality 
class should be made. Note, however, that the link between 
the centrality class and the impact parameter involves the total 
inelastic cross section, and to a very good accuracy for not-too- 
peripheral collisions c ~ 7rft 2 /<ri ne i. [631 - The total inelastic cross 
section is not measured, instead is obtained from Glauber simu- 
lations, thus carries model uncertainty 
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FIG. 7: The pionic HBT radii plotted as functions of the 
average transverse momentum of the pair compared to the 
STAR data |59j] at the centrality 0-5%. The calculation uses 
the two-particle method and includes the Coulomb effects. 
The values of the model parameters are the same as in Fig. \E\ 



The in-plane and out-of-plane freeze-out curves are de- 
fined as the intersections of the freeze-out hypersurface 
with the planes y = and x = 0. They are obtained 
with the initial central temperature Tj = 320 MeV and 
the final (freeze-out) temperature Tf = 150 MeV and are 
shown in Fig. [H The two freeze-out curves practically 
overlap, indicating that the expansion of the system is 
almost azimuthally symmetric in the transverse plane. 
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FIG. 8: (Color online) The freeze-out curves for 
Ti = 320 MeV, Tf = 150 MeV, and b = 7 fm. The solid line 
describes the in-plane profile, while the dashed line describes 
the out-of-plane profile. 



This effect is certainly expected for the almost central 
collisions, where the impact parameter is very small. We 
note that the shape of the isotherms is consistent with 
the result presented in Fig. 2 of Ref. [12], where the same 
decoupling temperature of 150 MeV was considered. 

In Fig. [6] we present our results for the hadron 
transverse-momentum spectra with the same values of 
the parameters. The dots show the PHENIX data [6(J 
for positive pions, positive kaons, and protons, while the 
solid lines show the results obtained from our hydrody- 
namic code connected to THERMINAT0R. Our model de- 
scribes the data properly up to the transverse-momentum 
values of about 1.5 GeV. For larger values of pt the model 
underpredicts the data. This effect may be explained by 
the presence of the semi-hard processes, not included in 
our approach. 

In our model, the values of T and b control the over- 
all normalization. On the other hand, the value of the 
freeze-out temperature Tf determines inter alia the rel- 
ative normalization and the slopes of the spectra. We 
note that the correct slope for the pions and kaons is 
recovered at a relatively high value of Tf, This is pos- 
sible, since the spectra of the observed hadrons contain 
the contributions from all well established hadronic reso- 
nances. This single-freeze-out picture (3ll If33l. I(34| . where 
the same temperature is used to describe the ratios of 
hadronic abundances and the spectra, was first tested for 
hadronic spectra in Ref. [63| , see also [H, [66| ■ We note 
that the value Tf = 150 MeV agrees with the recent re- 
sults obtained in the framework of the single-freeze-out 
model in Ref. - In the calculation of the proton spec- 
tra, in order to be consistent with the PHENIX experi- 
mental procedure, we removed the protons coming from 
the A decays. 

In Fig. the model results and the STAR data 
for the HBT radii are presented. Again, for the same 
values of the parameters, a quite reasonable agreement is 
found. Discrepancies at the level of 10-15% are observed 
in the behavior of the i? s ido, which is too small, and i?i on g, 
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FIG. 9: (Color online) Transverse-momentum spectra of tt + , 
K + , and protons. The PHENIX experimental results [601 ] for 
Au+Au collisions at ^/sjvjv = 200 GeV and the centrality 
class 20-30% (points) are compared to the model calculations 
(solid lines). The values of the model parameters are the same 
as in Fig. 




FIG. 10: (Color online) The elliptic flow coefficient vi- The 
values measured by PHENIX 0| at ^/FKw = 200 GeV and 
the centrality class 20-40% are indicated by the upper (pions 
+ kaons) and lower (protons) points, with the horizontal bars 
indicating the pr bin. The corresponding model calculations 
are indicated by the solid lines, with the bands displaying the 
statistical error of the Monte-Carlo method. The parameters 
are the same as in Fig. 
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FIG. 11: The pionic HBT radii plotted as functions of the 
average transverse momentum of the pair compared to the 
STAR data [H at centrality 20-30%. As in the case of the 
central collisions, the calculation uses the two-particle method 
and includes the Coulomb effects. The values of the model 
parameters are the same as in Fig. 



which is too large, probably due to the assumption of 
strict boost invariance. The ratio -R ou t/-Rsidc — 1-25 is 
larger than one, which is a typical discrepancy of hydro- 
dynamic studies. Nevertheless, when compared to other 
hydrodynamic calculations, our ratio R ou t/Rside is sig- 
nificantly closer to the experimental value. The ratio is 
almost constant as a function of k T , contrary to the de- 
crease observed in the data. 



B. Non-central collisions 

Next, we consider the centrality classes 20 - 30% and 
20 - 40%. The data are compared with the model re- 
sults obtained with the impact parameter 6 = 7 fm. The 
values of the initial central temperature and the final 
temperature are the same as in the case of the central 
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FIG. 12: (Color online) The freeze-out curves for central LHC 
collisions; b = 1 fm, T, = 450 MeV, and T f = 150 MeV. Sim- 
ilarly to the central RHIC collisions, the two curves overlap 
indicating that the system at freeze-out is symmetric in the 
transverse plane. 



collisions, i.e., T t = 320 MeV and T f = 150 MeV. The 
two freeze-out curves are shown in Fig. [H We observe 
that the out-of-plane profile is wider than the in-plane 
profile. This difference indicates that the system is elon- 
gated along the y axis at the moment of freeze-out. This 
feature is in qualitative agreement with the HBT mea- 
surements of the azimuthal dependence of -Rside- The 
comparison of the experimental and model transverse- 
momentum spectra is presented in Fig. [9l In this case 
we also find a good agreement between the data and the 
model results up to the transverse-momentum reaching 
1.7 GeV. 

In Fig. [inithe model results and the PHENIX data [ED] 
on the elliptic flow coefficient V2(pr) are compared. We 
observe that the V2 of pions+kaons is about 10% below 
the data. Taking into account the uncertainly in the ini- 
tial eccentricity, discussed at the end of Sect. Hill which is 
at the level of 10-20%, the obtained agreement is reason- 
able. On the other hand, the model predictions for v-i of 
protons is too large. The discrepancy is probably caused 
by the final-state elastic interactions, not included in our 
approach [68j |. 

The HBT radii for non-central RHIC collisions are 
shown in Fig. [TTJ Similarly to the central collisions, 
we observe that i? S idc is slightly too small and i?i ong is 
slightly too large. Still, the ratio i? ut /Rside is very close 
to the data. Comparing our values of i? s idc with other 
hydrodynamic calculations we find that our values are 
larger. This effect is caused by the halo of decaying res- 
onances which increases the system size by about 1 fm, 
see Ref. [H. 



VII. PREDICTION FOR THE LARGE HADRON 
COLLIDER 

We expect that increasing the initial beam energy re- 
sults essentially in a higher initial temperature Tj used 



FIG. 13: (Color online) The model results for the transverse- 
momentum spectra of 7r + , K + , and protons. The values of 
the model parameters are the same as in Fig. [12] 



as the input for the hydrodynamic calculation (in the 
midrapidity region studied here). Therefore, to make 
predictions for the collisions at the LHC energies we use 
a set of values for Ti which are higher than that used at 
RHIC, namely, % = 400, 450, and 500 MeV. It is not 
clear which value of Ti will be "selected" in LHC. Esti- 
mates based on extrapolations [3] suggest an increase 
of multiplicity by a factor of 2 compared to the highest 
RHIC energies, which would favor T around 400 MeV. 
However, to investigate a broader range of possibilities, 
we consider much higher temperatures as well. As alredy 
mentioned, in the case of LHC we also use a larger value 
of the nucleon-nucleon cross section in the definition of 
the initial conditions (a = 66 mb) and different values 
of the chemical potentials in the hadronization made by 
THERM I NAT OR (fi B = 0.8 and n s = /Xj 3 =0). 

In Tables [J and HU we list our results for the following 
quantitites: the total 7r + multiplicity dN/dy, the inverse- 
slope parameter for pions A, the pion+kaon elliptic flow 
V2 at pr = 1 GeV, and the three HBT radii calculated 
at the pion pair momentum kr = 300 MeV. The inverse- 
slope parameter given in Tables HI and Ull is obtained from 
the formula 



A = - 



dpi 



In 



dN„ 



2-KpxdpTdy 



(26) 



The tables show several expected qualitative features. 
Obviously, as the initial temperature increases, the mul- 
tiplicity grows. This is due a larger initial entropy, which 
causes a larger size of the freeze-out hypersurface. We 
find that the following parameterizations work very well 
for the multiplicity of ir + at LHC for the two considered 
centrality cases: 



dN 

dy 
dN 

dy 



12600(T 4 /GeV) 



.3.4 



6870(T l /GeV) 



3.4 



(b = 1 fm) (27) 
(6 = 7 fm). 



The power 3.4 works remarkably well. This behavior re- 
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FIG. 14: The model results for the pionic HBT radii. The 
calculation uses the two-particle method and includes the 
Coulomb effects. The values of the model parameters are 
the same as in Fig. [12] 



fleets the dependence of the initial entropy on T as shown 
in top left panel of Fig. [21 where for the relevant temper- 
ature range of 300 — 500 MeV we have the approximate 
scaling s/T 3 ~ T 0A . 

Similarly, for the slopes in the studied domain we have 

A(l GeV) = 0.629T, + 0.034 GeV, (6 = 1 fin) (28) 
A(l GeV) = 0.6291; + 0.044 GeV, (6=7 fm). 

The interesting feature is a constant (i.e. independent 



T t [MeV] 


dN 
dy 


A [MeV] 




idc [fm] 


-Rout [fm] 


-Rlong [fm] 


320 


250 


235 


4.3 


5.4 


6.1 


400 


558 


286 




5.3 


6.0 


7.6 


450 


837 


318 




5.8 


6.2 


8.0 


500 


1193 


348 




6.3 


6.5 


8.6 



TABLE I: Central collisions at RHIC (the second row) and 
LHC (the three lower rows): A set of our results obtained 
for b = 1 fm and four different values of the initial tempera- 
ture: Ti = 320, 400, 450, and 500 MeV. The columns contain 
the following information: dN/dy - the total pion multiplic- 
ity (positive pions only), A - the inverse-slope parameter for 
positive pions at pr = 1 GeV, -Rside, out, long - the three HBT 
radii calculated at the average momentum kr = 300 MeV. 



of Ti) shift by about 10 MeV from b = 1 to b = 7. The 
HBT radii increase rather moderately with Tj, as can be 
seen from the tables. 

In the following Sections we show more details, dis- 
cussing the spectra, t>2, and the HBT radii obtained for 
the initial temperature T, = 450 MeV. 



A. Central collisions 

In Fig. [I2]we show the freeze-out curves obtained from 
our hydrodynamic code with Ti = 450 MeV. Compar- 
ing to the corresponding central RHIC collisions with 
Ti = 320 MeV from Fig. we observe that the differ- 
ence in the initial temperature results in the longer time 
of the hydrodynamic expansion and a larger transverse 
size (both increase by about 3 fm). On the other hand, 
similarly to the RHIC results, we find that the two freeze- 
out profiles overlap, hence the system at freeze-out is, as 
expected, azimuthally symmetric in the transverse plane. 
We also note that the shape of the isotherms is consistent 
with the result presented in Fig. 4 of Ref. fl2l |. 

In Fig. [13] we give the model transverse-momentum 
spectra of hadrons. Compared to the RHIC results from 
Fig. [6] we find much larger multiplicities of the produced 
hadrons and smaller slopes of the spectra, indicating the 
larger transverse flow that is caused by the larger initial 
temperature. 



Ti [MeV] 


dN 
dy 


A [MeV] 


TT + K 


-Rsidc [fm] 


-Rout [fm] 


-Rlong [fm] 


320 


131 


245 


0.11 


3.8 


4.6 


5.6 


400 


303 


294 


0.11 


4.7 


5.2 


6.7 


450 


455 


327 


0.11 


5.1 


5.2 


7.1 


500 


651 


358 


0.11 


5.5 


5.3 


7.5 



TABLE II: Non-central collisions at RHIC (the second row) 
and LHC (the three lower rows) (b = 7 fm): The same quan- 
tities shown as in Table [J with the additional information on 
the pion elliptic flow «2 at pr = 1 GeV. 
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FIG. 15: (Color online) The freeze-out curves for pe- 
ripheral collisions at LHC, 6 = 7 fm, Tj = 450 MeV, and 
Tj = 150 MeV. The solid (dashed) line shows the in-plane 
(out-of-plane) profile. 




FIG. 17: (Color online) The elliptic flow coefficient v- 2 . The 
parameters are the same as in Fig. 1151 




Pt [GeV] 

FIG. 16: (Color online) The model results for the transverse- 
momentum spectra of 7r + , K + , and protons. The values of 
the model parameters are the same as in Fig. [15] 



such effects in data. 

In Fig. [16] we show the model transverse-momentum 
spectra of hadrons for the same values of the parame- 
ters. Compared to the RHIC case with T = 320 MeV 
and the same values of Tf and b, we find flatter spectra 
with higher multiplicity. In Fig. [17] we show our results 
for the elliptic-flow coefficient. The stronger transverse 
flow generated in this case induces larger splitting be- 
tween the pion+kaon vi and the proton i>2 , with the val- 
ues of the pion+kaon elliptic flow very similar to those 
found in the case T = 320 MeV. This result indicates 
the saturation of the elliptic flow of light particles for a 
given initial space asymmetry. On the other hand, the 
proton elliptic flow is significantly reduced. This observa- 
tion is consistent with the findings of Kestin and Heinz 
discussed in Ref. Finally, in Fig. [18] we show our 

model calculations of the HBT radii. We note that the 
ratio i?out/-R s idc is very close to one. 



Our model calculations of the HBT radii are shown 
in Fig. [TH The increase of the central temperature 
from T = 320 MeV to T = 450 MeV makes all the radii 
moderately larger. The ratio -R ut/-Rsidc decreases by 
about 10% which is an effect of the larger transverse flow 
caused, in turn, by the larger initial temperature. 



B. Non-central collisions 

In this Section we present our results describing the 
peripheral collisions with T = 450 MeV, T f = 150 MeV, 
and b = 7 fm. In Fig. [15] we show the freeze-out curves. 
One can observe that the system is initially elongated 
along the y axis, but in the end of the evolution it be- 
comes elongated along the x axis. This behavior is indi- 
cated by the crossing of the freeze-out curves. The change 
of shape is caused by the strong flow which transforms 
the initial "almond" into a "pumpkin" [fjjj]. Azimuthally 
sensitive HBT probes such non-trivial behaviour and can 
be used as a precise confirmation test for the existence of 



VIII. CONCLUSIONS 

Summarizing the results, we note that the increase of 
the initial temperature T of the hydrodynamic evolution 
yields a rather smooth change of the basic soft-physics 
observables, summarized in formulas ([27]l - 1(28)1 . The 
observed multiplicities, spectra, and HBT radii reflect the 
increased values of entropy and collective flow. One may 
of course wonder on the credibility of the predictions, 
which assume that the basic picture at LHC resembles 
RHIC, with the initial conditions cranked up. The LHC 
results in the coming years will answer that question. Of 
course, one may undertake extrapolations only when the 
current experimental situation is described well with the 
given model, which is not an easy task at RHIC, where 
no approach describes uniformly and accurately all the 
available data, including the femtoscopy. 

In this paper we used the hydrodynamic approach 
linked to the statistical hadronization to globally describe 
the midrapidity hadron production at the highest RHIC 
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energies and to make predictions for the future experi- 
ments at LHC. We have succeeded in the overall quite 
correct description of the soft hadron production at the 
highest RHIC energies, making fits of the pion spectra, 
i>2 and HBT radii. For all observables we have reached 
agreement at the level of 10-15%, which should be consid- 
ered satisfactory baring in mind the approximations used. 
We have provided uniform fits to all basic soft-physics 
information, including the pr-spectra, the elliptic flow, 
and femtoscopy. We have made a best effort with the 
presently available tools, that is the perfect-fluid hydro- 
dynamics with the state-of-the-art equation of state, the 
hadronization including all resonances, as implemented 
in THERMINAT0R, and the two-particle method with the 
Coulomb correction for evaluating the HBT radii. Hav- 
ing fitted RHIC at its highest energy, we have then made 
predictions for LHC. This was achieved by increasing the 
initial temperature which is, in our opinion, the main 
physical effect expected at larger beam energies. Pos- 
sible improvements of the present approach include the 
incorporation of the elastic rescattering in the final state, 
putting viscosity into hydrodynamics, as well as provid- 
ing more accurate initial condition, including the initial 
evolution of partons and the fluctuations of the axes 
of the second harmonic of the shape generated by the 
Glauber Monte Carlo calculations. 



Finally, we note that THERMINAT0R may readily be used 
as a Monte Carlo generator in the energy domain of LHC 
also for such purposes as the detector modelling and test- 
ing. 



FIG. 18: The pionic HBT radii. The values of the model pa- 
rameters are the same as in Fig. [15] 



We are grateful to Piotr Bozek for useful discussions. 
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